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I.  INTFOCHCTION 


are  using  automated  Elation  to  aerate  approve 

solutions  to  the  prognostic  equations  of  meteorology,  these  equations  are 
treated  in  the  form  that  would  arise  by  means  of  medal  analysis  and  truncation 
consequently  the  equations  take  the  form  of  coupled  non-linear  first-order 
ordinary  differential  actions,  the  ntmter  of  such  options  may  be  very 
large  if  many  modes  are  included  in  tto  analysis.  *  have  been  principally 
concerned  with  the  foliation  o,  a  general  method  in  the  work  reported  here. 

this  report  officially  covers  the  period  from  Dec.  20,  1971  through 

►tech  31,  1973.  Some  of  the  work  reported  here  reached  Its  natural  con¬ 
clusion  in  the  months  follodng  a„d  was  inflUimcea  ^  ^  ^ 

a  talk  on  these  matters  by  the  auttor  at  tte  Pand  Corporation  in  tfcrch  1973 

at  an  APPA  meeting  on  climatology.  In  the  ccxirse  of  this  first  fifteen 

months  we  have  made  seme  false  starts  as  we  attempted  to  orient  ourselves 

in  this  work.  We  do  not  report  on  these  here.  This  is  the  sole  written 

report  on  this  work;  work  along  these  lines  continues  ruder  APPA  Grant  . 
DA-ARO-D-31-124 . 

IWs  work  has  been  performed  in  conjunction  with  an  extensive  program 
to  investigate  climate  prediction  and  modification  and  consequently  we  are 
interested  in  the  long  term  behavior  of  the  atmosphere  under  the  influence 
of  the  driving  'force'  of  the  sun.  Hie  existence  of  the  yearly  cycle  of 
seasons  would  indicate  that  the  atmosphere  is  acting  as  a  steady  state 
system  under  the  influence  of  the  yearly  heating  cycle.  We  might  than 
envisage  that  climate  modification  will  occur  if  the  heating  effect  on  the 
atmosphere  shows  a  long  term  change;  naturally  if  *  ala3  ^  into  acoount 
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the  hydrological  cycle  other  forms  of  climate  modification  could  be  anticipated. 


The  principal  accomplishments  durincr  the  period  of  this  grant  have  been 
1*  Programs  useful  for  the  analytic  solution  of  coupled  mode  eouations 
were  developed  and  tested.  These  programs  were  written  in  PL1/F0FMAC.  ^ 
Sane  have  been  tra.  elated  into  MACSYMA.  ^ 

2.  Several  analytic  methods  were  investigated  for  their  possible  bear¬ 
ing  upon  the  problem. 


During  the  tern  of  this  grant  the  principal  investigator  was  the  only 
senior  participant  in  this  work;  two  graduate  students,  Michael  Polcari 
(Feb.  '72  -  Jan  *73)  and  Theodore  J.  Povcraft  (Feb.  '73  -  )  were  supported  by 
this  grant. 

The  mathematical  problem  that  we  face  is  the  solution  of  coupled  non¬ 
linear  equations  containing  explicit  driving  terms  whose  diurnal  frequency 
is  large  compared  to  that  of  natural  oscillations  of  the  driven  system.  The 
envelope  of  the  amplitude  of  the  driving  terms  have  a  slow  time  seasonal 
oscillation  and  perhaps  a  longer  time  scale  variation  because  of  variations 
in  the  sun's  output.  In  particular  the  natural  modes,  i.e.  the  cyclones  and 
anticyclones  of  the  middle  latitude  disturbances,  have  frequencies  corres¬ 
ponding  to  the  period  of  about  a  week.  Inasmuch  as  the  large  scale  turbulence 
accounts  for  the  bulk  of  the  energy  transport  to  the  poles  we  cannot  ignore 
the  presence  of  the  high  frequency  modes  in  the  solution  of  the  differencial 
equations.  As  is  well  known  it  is  possible  because  of  the  non-linear  nature 
of  the  physics  and  the  equations  that  high  frequency  terms  can  have  a  low 
frequency  effect.  Before  we  can  attempt  to  solve  the  steady  state  oroblem  we 
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must  in  sane  wav  account  for  effects  of  the  hiqh  frequency  modes  -  that  is, 
the  baroclinic  instabilities .  Barring  our  ability  to  solve  the  equations 
exactly  the  average  effects  of  the  high  frequency  modes  must  be  calculated 
by  some  approximation  technique.  One  basic  difficulty  is  that  the  satura¬ 
tion  level  of  the  excitation  of  these  natural  modes  is  not  given.  In  the 
approach  taken  here  we  will  replace  the  non-linear  effects  of  the  oscillatinq 
inodes  by  constant  and  linear  terms  in  the  differential  equation.  This 
approach  is  closely  related  to  the  introduction  of  'eddy'  viscosity  and  of 
course  the  physical  ideas  are  similar. 

There  are  two  interconnected  problems  that  are  of  interest: 

a)  study  of  instabilities  to  determine  the  level  of  saturation  • 

b)  study  of  steady  state  behavior  under  driving  forces. 

The  values  of  the  dynamic  variables  about  which  we  linearize  the  equations 
is  not  chosen  arbitrarily  to  be  zero  but  rather  by  inspection  of  the  results 
of  nimerical  integration  of  the  equations  or  by  inspection  of  the  equations 
themselves . 

The  analytical  vork  to  be  done  by  ocmputer  is  hedged  bv  a  number  of 
inherent  restrictions  beyond  the  obvious  ones  that  the  results  equal  or 
approximate  the  correct  solution.  Among  these  restrictions  are  the  following: 

a)  It  must  be  possible  to  specify  the  steps  to  be  performed  in 
algorithmic  fashion  so  that  the  procedure  can  be  mechanized. 

b)  The  number  of  terms  that  are  generated  should  not  be  so  large  that 
no  conclusion  could  be  drawn  from  an  inspection  of  the  analvtical  results. 

c)  The  calculation  should  not  be  a  numerical  evaluation  in  disguise. 

The  results  of  a  high  order  perturbation- theoretic  calculation,  for  example, 
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involvinq  manv  coupled  inodes  will  not  satisfy  conditions  (b)  and  (c)  .  We 
expect  that  a  combination  of  numerical  and  analytic  methods  mav  prove  rmst 
useful;  the  numerical  evaluation  of  various  analytic  terms  will  serve  to 

delineate  those  parts  of  the  multi-termed  analytic  answer  which  are  most 
important. 

The  method  we  employ  is  a  variant  of  the  well-known  Bogoliubov-Krylov- 

Mitropolsky  (BKM)  scheme. (3)  This  method  is  used  to  generate  the  so-called 

slaw-time-scale  equations  for  the  amplitude  and  the  frequency  of  a  non-linear 

oscillation.  The  starting  point  of  the  entire  procedure  that  we  often  follow 

is  numerical  calculations  from  which  we  can  determine  a  reasonable  beginning 
approximation . 

In  tie  next  section  we  shall  discuss  the  mathematics  of  the  method  we 
employ  and  its  implementation  to  date.  In  Section  III  we  discuss  a  model 
for  the  saturation  of  the  baroclinic  instability  which  illustrates  some  of 
the  ideas  worked  out  in  the  previous  section.  In  Section  IV  we  discuss  the 
application  to  a  system  of  three  couoled  modes  whose  free  oscillations  are 
described  by  elliptic  functions  and  some  preliminary  results  on  the  theory 
of  rotating  basin  phenomena  as  formulated  by  E.  N.  Iorenz. (3) 
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II.  Mathematics  and  Implementation 


A.  Mathematics 


Our  basic  goal  is  to  use  the  information  obtained  by  linearization  of 
the  equations  and  by  fjerturbative  expansion  of  the  solutions  to  determine  a 
replacement  for  the  non-linear  terms  in  the  differential  equations  in  the 
form  of  constants  and  of  terms  linear  in  the  mode  amplitudes.  This  is  the 
gist  of  the  BKM  method  of  equivalent  linearization.  Once  this  step  has 
been  made  the  resulting  equations  should  be  quite  easv  to  solve  by  use  of 
standard  techniques  used  in  the  analysis  of  linear  svs terns .  The  equations 
that  we  consider  are  of_the  form 


dB. 

l 

dt 


rjkiBkVDi 


(t) 


(II-l) 


where  the  B.  are  the  mode  amplitudes,  the  M.  .  and  tl^  r.  .  are 
1  ID  ilk 

constants,  and  the  are  the  driving  terms  which  are  explicit  functions 
of  time. 


In  general  we  will  be  looking  for  periodic  solutions  of  the  coupled 
mode  equations  III-l  by  perturbational  analvsis.  Often  the  results  of 
linearization  will  lead  to  growing  or  decaying  solutions  rather  than  oscil¬ 
lating  ones.  In  such  cases  we  mav  have  to  chanqe  the  values  of  the  B 

1 

about  which  one  linearizes.  (For  an  alternate  approach  see  the  next  section, 
(Eq.  ITI-6)  ) .  For  example  let  us  consider  one  of  the  most  famous  nonlinear 
coupled  mode  equations, the  so  called  predator-prev  equations.  ^  Let  x 
and  y  denote  the  values  of  two  populations, namely  a  prey  population  x 
and  a  predator  population  y  .  The  prey  x  will  increase  exponentially 


-5- 


except  that  the  predators  eat  it  up  at  a  rate  proportional  to  the  produce  of 
the  two  populations. 

In  symbols 

dx 

vr  =  ax  -  Byx 


Similarly  the  predators  would  decrease  exponentially  if  there  were  no  prey 
but  increases  in  proportion  to  the  number  of  prey  encountered.  So 


d y 

dt 


-  yy  +  6  xy 


we  can  make  the  identifications 


>"n  “  - 


*72  =  -  Y 


=  -  B 


If  we  were  to  linearize  these  equations  about  the  point  of  equilibrium 
x  =  y  =  0  ,  we  would  find  that  the  prey  increases  as  exp-At  and  the  predators 
decrease  as  exp-Bt.  Carrying  out  perturbation  theory  will  not  ]ead  us  away 
from  real  exponentials  whereas  the  actual  populations  oscillate.  We  can 
obtain  this  behavior  by  linearizing  about  a  second  point  of  equilibrium 
yo  =  o/B  ,  XQ  =  y/&  .  Writing  Y  -=  By  -  a  ,  X  =  6x  -  y  , 
we  get 


dX 

dt 


=  -  yY  -  XY 


dY 

dt 


=  aY  +  XY 
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If  one  linearizes  this  pair  of  equations  about  the  point  X  =  Y  =  0  the 
resulting  set  has  an  oscillatory  solution  at  the  frequency  2irf  = 


In  any  event,  based  upon  the  results  of  the  numerical  integration  arri/or 
seme  analytical  study  of  the  set  of  Equations  III-l  we  introduce  node  ampli¬ 
tudes  u.  by  =  ui  +  Bio  we  then  have: 

du. 

1 

dt 

where: 


I 


M.  .  u.  +  e  , 

13  3  jTk  rjk 


I 


U. 


\  +  Di 


(II-2) 


M.  .1 
ID 

=  M.  .  + 
ID 

?  (r. .. 1  b  + 
iDk  k,o 

rkfj 

D.1 

* 

-  D.  + 

Zh..B.  +  l 

r . 

1 

l  4 

j  13  30  jk 

jk 

ho 


We  have  introduced  a  purely  formal  ordering  parameter  e  which  is 
considered  small  for  the  purposes  of  carrying  out  a  perturbation  expansion 
but  which  in  actuality  is  equal  to  1. 


In  order  to  make  the  analysis  seem  less  abstract  we  carry  out  the  steps 
indicated  symbolically  on  a  specific  example,  viz.  a  driven  harmonic  oscil¬ 
lator  on  which  a  fictional  force  of  the  form  ev|v|acts  (v  is  the  velocity) . 
Ihe  equations  governing  the  notion  of  this  oscillator  are 


dv  2 
dt  =  "  u  x 


ev|v|  +  D  sinftt 
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— — 


Ihe  next  step  is  to  find  the  eigenmodes  of  the  linear  part  of  (II-2) 

and  then  make  a  linear  transformation  so  that  the  coupled  node  equations 
read: 


d6. 


1  "  "i 

jr-  =  y.$.  +  D.  +  er..  0.6, 
at  111  jk  3  k 


(II**3) 


in  which  the  0's  ,  the  D  '  s  and  the  r  '  s  are  linear  combinations  of  the 

I  I 

p  's,  the  D  's  and  the  r's  .  In  the  case  of  our  example  the  equations 


are 


dt  {x  +  “  iaj  (x  +  VM  +  r-  sinfit 


D 

iw 


and  its  octiiplex  conjugate. 


We  denr’-e  the  solution  of: 


d3. 

ar "  wi3i +  Di 


(II-4) 


by  I£W  (I)  where: 


LOW(I)  =  0i(o)e 


p.t  t  y. (t  -  t')  „ 

+  /  e  D_.  (t')dt' 


(II-5) 


once  again  our  exanple  would  give 


v 


■  .  iut  t 

“w  -  +  5—  /  D  sin£2t'e”lut  dt' 


From  this  we  write  the  next  tern  in  the  formal  perturbation  expansion  of  the 
solution  of  (II— 4 )  in  the  next  order  as: 


-8- 


(2) 


t  y. (t  -  t ' )  V 

=  /  e  Z.  r1 


jk  3^ 


I£W  ( J)»ICW  (K)  dt ' 


(H-6) 


Inserting  the  lowest  order  results  into  this  equation  we  have: 


(2)  V 

*i  «/  exp[y.(t  -  t')  ]  Z.  exp[(v.  +  y.  )t']  . 

-!  *  J  K 


jfk 


(II-7) 


Fjk  Bj(o)  +  fQ  e  j  Dj(x)  dT  *  Bk(0)  +/  e  Uk  D^(T)dT 


Next  we  can  substitute  the  expansion 


V  •  —  u. 

1  p  lpe 

into  (II-7)ard  then  carry  out  the  integrations  straight- forwardLy.But  if  we  were 
to  do  so  a  classic  difficulty  might  be  encountered; namely  that  the  frequency 
of  one  or  more  of  the  terns  on  the  right-hand  side  of  (II-6)  is  equal  to  y  . 
Such  a  term  is  called  secular  and  when  integrated  leads  to  a  term  that  goes 
linearly  with  time  and  hence  is  unbounded. 

We  discuss  this  situation  in  sene  detail  in  Appendix  A  and  show  how  such 
secular  terms  can  be  treated  by  appropriately  shifting  the  eigenvalues. 

After  this  has  been  done  one  can  carry  out  the  integrations  indicated  in  (II-7)  . 

Returning  to  the  example  of  the  nonlinearly  damped  oscillator  that  was 
introduced  above  we  shall  apply  the  procedures  (the  BKM  method)  given  in 
Appendix  A.  First,  consider  the  case  of  the  undriven  oscillator.  We  have 
then,  with  A  =  x  +  v/(ioo) 


=  iuiA  -  v  |  v  I 
at  loo  1  1 
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Hie  lowest  order  solution  is  then  A  =  e1*,  4,  =  ut  .  Assuring  that 
A  =  ae1^  +  eu (a,  \p)  +  ...  and  that 

dii 

—  a)  +  c  (a)  +  •  • 

•  da  „  .  . 

gj-  =  e  Cj^  (a)  +  . . . 

one  finds  that  Bj,  =  0,  •-  -  -  a2w  .  Here  we  have  used  the  BKM  method. 

Ihe  implied  change  in  the  decay  rate  of  the  lowest  order  solution  is  the 
device  used  to  counter  the  terms  proportional  to  eilJ;  would  arise  upon 
inserting  A  =  ae1^  and  its  complex  conjugate  into  the  term  -ev|v| . 

This  procedure  can  be  carried  one  step  further  by  means  of  the  so -ceil  led 
method  of  equivalent  linearization.  The  results  that  have  just  been  obtained 
indicate  that  the  amplitude  damps  at  a  rate  given  by 

1  da  4 

a  dt  T  zabi 

Let  a  denote  the  average  amplitude  corresponding  to  some  given  initial 
condition.  Then  we  write  the  equation  for  a  linear  oscillator  that  is  in 
some  sense  equivalent  to  the  nonlinear  one  as : 

dv  8  -  ,2 

gjr  +  j  eawv  +  u  x  =  0 

The  solution  of  this  equation  is,  to  first  order  in 

x  =  x(0)e  ea*"cos(wt  +  <fi)/cos  (<f>) 
with  4  to  be  determined  by  initial  conditions. 
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In  the  case  of  the  driven  oscillator  we  see  that  the  effective  damping 
rate  will  depend  on  th;  driving  forces  since  these  determine  the  average 
amplitude .  If  we  were  to  ignore  the  nonlinear  term  and  look  for  the  steady 
state  of  the  oscillation  under  the  driving  force  we  would  find  that 


D  sinftt 


Taking  a  =  D/|^2  _  j  we  could  write  the  complete  equation  for  the  oscil¬ 
lator  as 

dv  ,  8  -  ,2  _  .  ^  ||,8  - 

gt  +  J  ea-^  v  +  w  x  =  D  sinfttrev|  v|  +  -j  eaP-  v 

TT  TT 

The  method  of  equivalent  linearization  chooses  the  coefficient  of  v  on  the 
LHS  of  the  above  equation  just  so  that  the  effects  of  the  last  two  terms  on 
the  right  hand  side  of  the  equation  cancel  out  when  averaged  over  a  cycle. 
Then,  in  fact,  if  one  were  interested  only  in  the  response  at  the  frequency 
n  one  would  write 


dv  ,  ,8  eDP  x 

dt  '3  2  02 

0)  -  ft 


v  +  o)  x  =  D  sinftt 


The  principle  result  from  the  analysis  has  been  the  determination  of  an 
effective  linear  damping  coefficient,  one  that  is  dependent  on  the  amplitude 
and  the  frequency  of  the  driving  forces. 
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B.  Implementation 


Almost  all  of  the  analytical  processes  mentioned  in  Section  A  have  been 
implemented  by  programs  written  in  PL/1  FOFMAC.  This  includes  the  processes 
of  calculation  of  the  secular  determinant  of  the  linear  part  of  the  system, 
linearization,  diagonalization,  perturbation  expansion,  and  substitution  of 
the  second  order  results  tack  into  the  differential  equations.  In  essence 

only  the  last  steps  involved  in  eigenvalue  re-evaluation  have  not  been  com¬ 
pletely  implemented. 


The  listings  of  these  programs  is  given  in  Appendix  B.  The  program 
names  and  their  functions  are 

FIRST  shifts  the  points  of  linearization. 

SBOOND  Ffnds  eigenvalues  and  left-and-right 

eigenvectors . 


THIRD 


Carries  out  diagonalization  formation 
of  lowest  order  solution  and  parts  of 
perturbation  expansion. 


Of  these  three  programs  FIRST  and  THIRD  have  been  translated  into  MACSYMA. 

It  is  probably  better  to  perform  the  function  of  SECOND  by  a  numerical 
program  since  the  eigenvalues  cannot  be  found  analytically  very  readily  ever 
for  a  system  with  five  modes.  The  programs  run  as  listed;  they  are  not 
complete  in  the  sense  that  our  ideas  are  still  developing. 


We  have  also  "boiler-plated"  a  program  together  which  numerically  inte¬ 
grates  ordinary  differential  equations,  plots  the  solutions,  and  also  Fourier 
analyzes  these  solutions.  We  have  used  this  program  to  investigate  the 
qualitative  nature  of  the  solutions  so  that  we  can  begin  the  analytical 
approximation  of  the  solutions. 
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Model  for  Saturation  of  the  Baroclinic  Instability. 


In  this  section  a  pair  of  coupled  nonlinear  equations  describing  the 
time  evolution  of  a  baroclinic  instability.  This  model  is  completely  solvable 
and  we  can  illustrate  how  our  perturbative  approach  weald  succeed  in  approxi¬ 
mating  the  correct  solution.  We  use  the  same  model  as  Saltzman  and  Tang^  do 
in  their  perturbative  calculation  of  this  instability  including  these  features: 


1)  a  two  level  model  on  a  beta  plane  of  finite  width  (a  beta- 
strip)  is  employed,  with  x  in  the  east-west  direction 
and  y  measured  north- south. 

2)  the  thermal  wind  relation  is  taken  to  be  fip  =  $ 

3)  we  look  for  a  wave  solution  of  the  form 


*1 

*3 

“2 


} 


v  e 


ikx 


sin  Zy 


which  vanishes  at  the  two  edges  of  the  beta-strip. 

4)  the  base  flow  has  a  thermal  wind  given  by  UT  =  (u^  -  U^)  /2 
and  a  mean  flew  given  by  UM  =  (L^  +  u3)/2  both  of  which  are 
in  the  east-west  direction. 


The  equations  governing  the  first  order  eddy  fields  (the  stream  function) 
can  be  written  as 


d 

dt 


=  ik 


-RM 

(2r-l)UT 


III-l 
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W^lPXe  “  ^1  +  ^3^/2  '  ~  "  ^3)/2  •  ^  and  Rj.  are  the  speeds 

of  the  two  Rossby  waves  which  are  the  solutions  of  these  equations  when 

^  =  0  •  quantity  r  is  given  by  [1  +  (k2  +  t^/X2]"1  where  X  is 
the  characteristic  reciprical  length  f/Salp  . 


*  = 


It  is  convenient  to  introduce  new  amplitudes  by  means  of  the  substitution 
+  %)/2]t  ^  then  dropping  the  primes.  One  gets 


1  1  R 

ikUT  dt  ^2r-l)  -  a 


*  V 

<tv 


(III- 2) 


where 


„  =  df/dy  r 

2  2  2 
i  +  k 


The  eigenvalues  of  the  matrix  on  the  right  hand  side  of  equation  III-2 
given  by 


X  = 


±  i/ 


(2r-l)  - 


2 

a 


Denoting  by  P  the  quantity  -  i>'  (2r_])  _  Q2  (so  that  ikP  UT  is  the 
positive  growth  rate  of  this  baroclinic  instability  in  the  limit  of  vanishing 
amplitude)  we  introduce  the  linear  combinations  <t>+  =  ^  (a  ±  P)  for 
which 


at  *±  -  *  ikuTp 

Na-;  the  resuits  of  the  S-T  calculation  in  second-order  show  that  there 
is  a  mean  zonal  flow  which  varies  as  sin(21y)  and  which  has  equal  and  opposite 
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values  on  the  first  and  third  levels.  If  we  denote  the  amplitude  of  this 
zonal  flow,  (call  it  the  2nd  harmonic  flow)by-2'f'  on  level  1  ard  -t-2'f  on  level 
3  then  we  find  that  (ignoring  ct  ccnpared  to  P  )  ’. 

d<|>+ 

ar~ =  (1  +  ^  *+ r)  (in-3) 

The  second  harmonic  flow  in  turn  is  fed  by  the  nonlinear  terms  in  the  vorticitv 
equation;  its  rate  of  increase  is  determined  by  the  product  iji  .  within 
the  context  of  the  perturbation  calculation  then  we  have 


df>+ 

ar  =  (v "  x 


dn-4) 


vhere  K  is  a  positive  constant  and  depends  on  the  other  parameters  (k,  jj, 
f,  etc.)  of  the  system,  v  =  kUT|p|  and  X  =  kU^/P2  .  m  this  last 
equation  we  should  include  terms  proportional  to  $  2  and  to  v+  <J>  ; 
however,  during  the  growth  of  the  instability  these  terms  will  be  relatively 

#  O 

small  in  comparison  with  4>+  . 

One  finds  the  solution  to  this  set  of  equations  to  be  given  by 


4>,  (0)  cosh  (/C  T) 

4+  =  — - — 

cosh  (»^  (T  -  t)) 

vtere  ,/c  i  v  and  e"2^  T  2;  XK(f>o2/8C  «  1 


(III-5) 

T  is  the  time  for  saturation. 
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is  large  we  can  write  the  solution  given  above  as 


Since 


+JC  T 
e 


vt 


x  'V  4*(0)S  _  x  l r\\  /_ 

*  -  -2vT  2vt  ~  *(0)  (e 

1  +  e  e 


vt 


-  e 


3vt  -2vT. 


+ 


•  ) 


(III-6) 


3  , 

for  t  <<  T  .  Ihe  seoond  term  'v  eJV  would  be  the  perturbation  result. 
Presumably  one  might  be  able  to  work  backwards  fran  perturbation  theory 
applied  to  a  growing  mode.  We  also  find  that 


Hf  =  £  (1  +  tanh  [v'c  (t  -  T)))  (III-7) 

and  f-  =  growth  rate  =  -  v  tanh  (t  -  T)  ]  (III-8) 

The  results  indicate  that  at  saturation,  i.e.  when  cty/dt  ->  0  one  has 


(III-9) 

max 

It  is  seldom  possible  to  find  an  exact  solution;  in  most  cases  one  must 
use  perturbation  theory.  Perturbation  of  a  growing  mode  yields  terms  that 
increase  faster  than  the  lowest  order  solution  (see  Equation  III-6  above  for 
an  interpretation  of  such  terms) .  In  terns  of  the  formalism  discussed  in 
the  previous  section  which  we  intend  to  use  for  perturbation  calculations 
one  rewrites  equations  (III-4)  in  terms  of  the  variable 


*+  “  *rnax/2 


+  ♦’ 


(III-10) 
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to  obtain 


IV.  Further  Examples 


We  illustrate  the  method  outlined  in  Section  II  by  means  of  a  relatively 
simple  but  non- trivial  problem.  Consider  three  modes  that  are  coupled  to¬ 
gether  and  which  are  driven  by  a  external  driver  in  the  following  fashion: 

dA 

gjT  =  Be  +  D  sin  (fit) ; 

dB 

dt  “  ”  ^  ;  (iv-i) 

|=-k2AB;  k2  <  1  ; 

In  the  absence  of  the  driver  the  solution  of  these  equations  (with 
A(0)  =  0,  B(0)  =  1,  C(0)  =  1)  would  be: 

A  =  sn(t,  k) 

B  =  cn(t;  k) 

C  =  dn  (t,  k) 

where  sn,  cn,  and  dn  are  the  usual  Jacobi  elliptic  functions.  We  consider 
the  driven  case  here.  As  a  first  step  we  integrated  this  set  of  equations 
numerically  using  the  initial  values:  A=0,B  =  1,  C  =  l.  The  numerical 
solutions  indicated  that  the  mode  amplitude  C  had  a  non-zero  average  value 
close  to  1  while  the  amplitudes  A  and  B  averaged  to  zero.  Therefore, 
we  linearized  the  equations  about  the  "point"  A  =  0,  B  •=  0,  C  =  1  using 
the  program  FIRST  (See  Section  II) .  The  matrix  of  the  linear  part  of  the 
new  equations  is 
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Next  we  solved  the  secular  determinant  and  formed  the  left  and  right 
eigenvectors  using  program  SECOND.  The  crank  was  turned  further  and  in 
Figure  2  we  show  seme  of  the  results  of  substituting  the  perturbation 
solution  back  into  the  differential  equation. 

Even  in  this  relatively  simple  example  we  find  that  very  many  terms 
oould  result  from  these  simple  calculations. 

The  BKM  method  permits  selection  of  those  terms  that  lead  to  changes 
in  the  growth  rate  and  frequency  of  the  oscillations  (waves) . 

In  order  to  test  "experimentally"  if  there  was  any  validity  to  the 
methods  we  have  proposed  in  Section  III  we  integrated  the  equations  IV- 1 
numerically. 


At  the  same  time  we  also  integrated  the  following  linearized  equations 
numberically. 

dt  =  1//>1  +  k2/4  B  +  D  sinDt 

dB  .  .  - k — 

dt  ”  "  1//vl  +  k2/4  A 

In  both  cases  D  =  .5,  k2  =  0.5,  and  0  =  0.2  . 

The  point  of  this  calculation  was  to  determine  if  the  solution  of 
jjggggiggd  equations  approximates  the  exact  solution.  (The  frequency  cor¬ 
rection,  given  by  the  reciprocal  square  root,  is  the  standard  first  order 
result  for  elliptic  functions) .  Figures  3  and  4  show  the  results  for  A  and 

B  ,  the  correspondence  is  sufficient  to  indicate  that  there  is  merit  in  the 
approach. 
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We  have  also  just  begun  an  application  of  the  methods  discussed  in 
previous  sections  to  the  modal  equations  derived  by  E.  N.  Lorenz  ^  to 
model  the  rotating  basin  experiments.  We  have  integrated  these  equations 
numerically  for  certain  values  of  his  frictional  and  hearing  parameters  as 
well  as  for  initial  conditions.  On  the  basis  of  these  numerical  results 
we  carried  out  seme  of  the  analytic  processes  described  in  Section  3  by 
means  of  the  programs  that  have  been  developed.  Figure  5  shows  the  results 
scxne  typical  numerical  results.  This  work  is  continuing. 
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APPendix  A-  glow  Time  Scale  Bjuatinnc; 

In  this  apperi-i:  we  shall  give  a  brief  description  of  the  Bogoliutov- 
Krylov-Mitropolsky  (BKM)  setae  for  deriving  slow  thus  scale  Rations 
Whose  solutions  are  "asymptotic"  to  the  solution  of  a  nonlinear  differential 
equation.  After  discussing  the  BKM  for  an  oscillator  with  one  degree  of 
freedom  in  the  conventional  manner  we  introduce  the  concepts  connected  with 
the  method  of  equivalent  linearisation.  The  use  of  omplex  notation  is 
then  introduced  as  Ons  is  more  natural  for  the  equations  we  consider. 
Finally  we  discuss  the  case  of  many  degrees  of  freedom 

Ihe  nonlinear  equations  that  we  look  at  first  are  of  the  sort  that  arise 
in  the  mechanical  oscillations  of  a  system,  viz 


m  o  +  ^  “  Ernf(x,  x) ;  x  r  • 
dt  dt  ' 


(A-l) 


Wherein  we  are  considering  a  system  with  one  degree  of  freedom  and  where  f 

is  a  (non-linear,  function  of  x  and  i  .  furthermore  is  a  parkier  which 

is  to  be  considered  small.  The  BKM  seta*  yields  solutions  in  powers  of  « 

which  are  as„totic  to  tie  solution  of  (A-l,  in  the  mathematical  sense,  i.e. 

a  finite  number  of  the  terns  in  the  expansion  well  approximates  the  true 
solution. 


becanes 


First  consider  a  simple  perturbation  solution  in  the  case  so  that  (A-l, 


d2x  2 
dt2  X  “ 


3  2  k 

-  ex  ,  <d  =  — 
o  m 


(A-2) 
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If  we  assume  that  we  can  write  the  solution  as 


>  n  (n) 
x  =  2L.  c  x  (t) 


(A-3) 


n=o 


we  obtain,  upon 
order  in  s 


inserting  Equation  (3)  into  (2;  and  equating  terms  of  equal 


+  „  V°>  =  o 

at2 

,.2  ° 

dt 


(A-4-a) 


(A-4-b) 


This  is  the  so-called  Bisson  method  and  it  leads  to  immediate  difficulties. 
The  general  solution  of  (A-4-a)  is 


X(0)  =  B  oos(wQt  +  4>) 

where  B  and  $  are  constants.  Now  then 


3  B  3 

x(°)3  =  cos(uot  +  4>)  +  cos (3a)Qt  +  34>) 


(1) 


This  means  that  the  right  hand  side  of  the  equation  for  x  has  a  tern 

(3B  }  4  •  oos(wot  +  »)  )  which  is  a  solution  of  homogeneous  equation.  As 
is  well  known  the  particular  integral  will  then  have  a  term  that  behaves  as 
t  sin  (cot)  .  This  is  the  so-called  secular  behavior  which  is  unbounded  in 


time. 


non 


The  physics  of  the  problem  illuminates  the  mathematical  solution.  Ihe 
-linear  restoring  forces  much  change  the  frequency  of  the  oscillation 
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because  the  potential  energy  is  not  simple  kx2/2  here  but 


rather  is 


PE  =  kx2/2  +  cmx4/ 4 


Consequently  the  solution  is  approximately  given  by 


B  oos (wQt  +  6ut  +  *)  =  B  cos(wQt  +  *)  -  «rt.  B  sin(u)  t  +  $) 


with  6u>  arising  from  the  non-linearity.  While  the  names  of  Lindstedt 
and  Poincare '  are  associated  with  methods  for  obtaining  period  solutions  of 
A-l  we  shall  describe  here  the  method  due  to  Bogoliubov,  Krylov  and  Mitropolsky; 
the  method  is  somewhat  heuristic  and  intuitive  but  paverful  and  direct. 


One  starts  from  the  fact  that  the  solution  to  A-l  with 


c  =  0  is 


a  cos  iJj  with  gr-  -  0  ,  (a  =  constant)  and 


We  recognize  tint  the  nonlinear  terms 


dt  “  wo  '  *  =  w0t  +  ♦  • 


—  —  nonlinear  terms  my  cause  a  change  in  the  amplitude, 
'  on  a  Slav  time  scale  and  in  the  frequency,  di|»/dt  .  We  write  then 


x  =  a  cos  ip  +  E  Uu;  (a,  $)  +  ... 

• 

x  =  -  ojQa  sin  ip  +  . . . 

dil;  (1) 

dt  =  wo  +  e  Cu;  (t)  +  ... 

da  _ (l)  .  . 

dt  “  e  A  (t)  +  . . . 


(A-5-a) 

(A-5-b) 


(A-5-c) 


(A-5-d) 


where  C  and  A  ;  are  to  be  determined  so  that  U(1)  will  contain 
secular  terms.  One  finds  then 
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(A-6-a) 


dip 


/2tt 

1 


dt  o  2ircoa  I  ^  (a  cos  ip,  -aw,  sin  \p)  oos  p 


di p 


^2ir 


da 

dt  2ttw  I  ^  (cos  \p,  —  a 
'o 


wQ  sin  ip)  sin  ip  dp 


Returning  to  the  example  we  started  from  in  which 


(A-6-b) 


f  (x,  x)  -  mx^  =  ma^  cos\ 


we  find 


dp 

dt  “  w 


2*1  f 

o  2irw  J 


cos  p  dp  = 


w  - 


ema~ 

2ttw 


da 

dt=° 


We  can  physically  interpret  the  results 
that  the  integral  in  (A-6-b)  is  the  work  done 
involves  the  so-called  reactive  power. 


given  by  A-6  by  recognizing 
per  cycle  while  that  in  (A-6-b) 
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B.  The  Method  of  Equivalent  Linearization 


There  is  another  way  of  interpreting  the  results  given  above  using  a 
somewhat  cruder  approximation.  Recall  that  we  started  by  looking  at 

d2x  2 

— 2  +  “  x  =  “  mcf (x,  x) 
dtz  ° 

The  lowest  order  BKM  solution  (the  eU(1)  term  can  be  ignored  here)  is: 


with 


x  =  a  cos  \p 


—  =  eA(1) 
dt  eA 


(a) 


We  approximate 


31  =„  +£C(1> 

at  o 


these  last  two  equations 


(a) 

as 


where  aQ  is  the  amplitude  for  t  =  o  .  One  observes  that  the  corresponding 

solution  for  x  could  have  been  obtained  from  the  solution  of  the  differential 
equation 


+  2Xx  +  —  x 
m 
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if  to  the  first  order  in  e 


we  had  chosen 


A 

k 

m 


2 

w  +  2a ) 
o  o 


eC'  (a) 


Wo  note  this  use  of  perturbation  theory  to  nodify  the 
by  linearization. 


eigenvalues  obtained 


C.  Complex  ftamulation. 


■me  previous  results  can  also  be  viewed  Iran  the  results 
matrix  equations,  lhe  second  order  differentia]  Equation  A-2 
in  matrix  foun  as  (use  P  =  x/m) 


Introducing  a  -  x  +  P/im^  and  its  oonplex  conjugate  a*  as 
we  obtain 


da 

or  =  in  a 
dt  o 


mu 


f  (c 


* 

a  ) 


o 


and  its  complex  conjugate.  Here  «  =  /kfa 

oo 

Writing  a  =  ae1^  +  ^  cn  6  {n)  (a,  *) 

1 


along  with 

oo 

at =  u0  +  X  E"  4'{n)(a) 

n**l 


da 

dt 


J  cn  A(n)  (a) 


one  obtains 


-  w  1_  - 
o  3i|]  c  >  - 


-  (A(n)  +  i?(n)a)  +  c-i* 


E*  -  ¥  (An-S  i-  +  Yn-S  5_)  B  <s> 
x  aa  9i y  p 
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on  first  order 
can  be  written 


netf  variables 


is  the  coefficient  of 


f" 


in  the  development  of  mef  (x, 


In  order  to  avoid  secularities 


one  must  have 


(n) 


if^a  = 


2tt 
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x)/iw  . 


■)  6(S) 


- - ,V~' 


D.  Several  Degrees  of  Freedom 


In  abstract  form  we  must  solve 


3i  C"  -A 

r  =  ZMij  Bj  +  Fi  (B)'  1  =  ^  ••• 


Itike  S  to  be  the  matrix  that  diagonalizes  M  ,  i.e.  S_1MS  =  A  where  A 


is  diagonal.  Take  U  =  S  .  B.  .  Then 

e  e]  ] 


dUk  .  > 
dt-=  1  Ak 


Again,  assuming  U  1  =  a  e11^  +  (a,  \p)  , 


5F=  xk  +  c^fa)  , 


da  T\  f  \ 

dt  =  c  Va) 


one  obtains 


*i  +  iTa  =  k 


for  the  eigenvalue  corrections  to  the  k-th  mode. 
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Appendix  B.  Program  Listings 


-31 


INPUT  TO  FORNAC  f rtZFFCCFSSCF 
FiFS  JP.OCEDUnE  C  £•  T I C.  N  2  ( M  A I  N )  ;  FORM  AC  C  ">  T 1 0  K  S  • 

/*  BASED  PARTLY  UPON  D^,CK  FOE  PERTU°H  t  o  [I-1 1  \  ^ 

/*  STARTED  JUNE  27  1972  *  / 

/*  COR  PEC  I  PEiiTUF.E  DECK  V 
CALL  IIIFUT; 

/*?.*/  CALL  LlNEAr.IZE; 

/*!,*/  CALL  FORM  DETERMINANT ; 

/*4.*/  CALL  0 U Q  P U T_ A L L ; 


/* 


SUBROUTINES 


V 


OPTSET (PRINT) ; 

INPUT : PROCEDURE ; 

/*!•*/  GET  FILE (CCNTi OL) DATA  (KCDES. 

NOTHING) ; 

LET  (FNC  (SIN  E)  =  -#I/2*  EXFON.  (#1*1  (1)  )  +  «  1/2*"  EXPON  .  (-#i*$  fin  • 
OPTSe?<;eXPND),;,=1/2"*:XPON*  ( #I  *' *  ( >  +  V2*  FXPON  .  (-»J*$1l))  )  ; 

/<2.>t'/  /*  ZERO  THE  ARRAYS  OF  COEFFICIENTS  * / 

/+AV  CALL  FOhEAC_L:,A5F  (’  D  (I)  «  ,  MODES,  1)  • 

/+HV  CALL  FOr.NAC_irAEE  (*M  (I, .7)  »  , MODES, 2)  ; 

/i  C*/  CALL  FCRUAC_F  RASE  ( '  G  AMU  A  ( I , »J  ,  K)  '  ,  MODES.  3)  * 

/*DV  CALL  FORM  AC_  ERASE  ('  80  ( I)  '  , MODES , 1 )  ;  ’ 

FOPMAC_EPAS  E: PROCEDURE (A, N,FANK)  ; 

DCL  A  C  H  A  F  ( *  )  /  ( *’•  /RANH  ,  N  i'  ( A )  )  FIXID  BINARY; 

DO  1  =  1  TO  4;  I r  1<=EANK  THEN  NT(I)=U;  FLSV  NT  (11=1*  run* 

ZEEO:  DO  1=1  10  N  *  (1,  ;  ChT  (I)  ;  DO  '.7  =  1  TO  K' '(?  cfr  (])' • 

DO  K*1  TO  1,1(3);  CEl'(K);  DO  1  =  1  t„  »-((.,  .  <U  n  .  '  1  ’ 

NET(A  =  0);  END  ZFiO;  A7CKIZE (I ; J ; K; L)  ;  END  FCPKACJJP  ASE  ; 

/*3.*/  /+  INPUT  I  HE  ARRAYS*/ 

/*A*/  CALL  EQUATIONS (« DRIVER' )  ; 

/+R*/  CALL  EQUATIONS ('KARFAY’ ) ; 

/*C*/  CALL  EQUATIONS (' GAMMAS ' ) ; 

FQU  AT  ION  3  ;  P  ROC  EDUP  E  ( A)  ;  DCL  A  CHAU(*)’ 

!  G0T0  F0F;  “•« 

^:^l^n^r^)COeY:  FP'”|V'LUF):  G0"J  wop; 


DCL  FOE  MAC 


EOS  R 


FOR  M  AC_  ER  Ab  E  ENTRY  (CHAR  (IS)  V  A  R  ,  F’  I X  E  D  BIN,  FIXED  R I  N ) 
EQUATIONS  ENTRY  (CHAR  (6)  VAR)  ;  '  >’ 

END  INPUT; 


LINEARIZE  ;  PF0CECUR1. ; 

IMODF.S:  DO  1=1  TC  MODES;  OEMI)- 
LET  (SUN  (1)  =0;  SDK  (2)  =0)  ; 
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J  MODES  :  DO  .1=1  VO  MOD  FT  ;  CET  (J)  ; 

LET  (SUM  (1)  =C)  ; 

DO  K—  1  TO  '1CDTS  *  CFT(F)  ; 

LET  (SUM  (1)=SUE  (i)  +  G\M'<A  (I,  J,K)  *  DO  (X)  +GAMVA  (I,K,J)  *BC  (K )  ; 
SUM  (2)  =SUM  (2)  ♦GAMMA  (I,  J,K)  *BO  (J)  '  BO  (K)  )  ; 

END; 

LET  (SUM  (2)  =SUK  (2)  +11  (I,  J)  ♦DO  ( J)  )  ; 

LET  (M  (1,0)  =M  (I,  J)  +SUM  (1)  )  ;  ATOMIZE  (SUM  (1)  )  ; 

END  J MODES  ; 

LET  (D  (I)  =D  (I)  +  S U h  (2)  )  ;  ATOMIZE  (SUM  (2)  )  ; 

END  IMODES; 

EOSK:  RETURN  ;  END  LINEARIZE; 


OUTPUT  ALL : PROCEDURE; 

/*AV  CALL  DISK  ('  D  (I)  1  ,  MODES,  1  DRIVE')  ; 

/*B*/  CALL  DISK  ('GAKiA  (IfJfK)  •  , MODES, 3, • GAM* ) I 

MOTEST= 1 1 • B  * 

/*C*/  CALL  Di SK  ( '  ROW  (I)  '  ,  MODES,  1 ,  '  ROWS'  )  ; 

RETURN; 

END  OUTPUT  ALL ; 


DCL  DISK  ENThY  (C.IAR(LO)  VAR, FIXED  BINARY, FIXED  BINARY, CHAR (20) VAE) 
DISK:  PR  OCE  DUKE  (A,  N,M,  B)  ;  OPEN  FIL17  (PUNCH)  OUTPUT  TITLE  (B); 

DCL  ( A  ,  D)  CHAR  (»)  ,  (N,i1,NT  (4))  FIXED  BINARY; 

DO  1  =  1  TO  4;  IF  I<=M  THEN  NT(I)=N;  ELSE  NT  (I)  =1  ;  END; 

CARDS: 

DO  1  =  1  TO  NT  ( 1 )  ;  CPS  (1)  ;  DO  1  TO  NT(2);  CE'?  (J)  ; 

DO  K=  1  TO  la  (3)  ;CLT  (K)  ;  DO  L=1  TJ  NT  (4)  ;  CET  (L)  ; 

IF  NOTES'"  THEN  GO  10  PUNCHER; 

IF  -«D  Ell  FMCG  ( 1  /  9999  1 97=  1  4  i  A ,  '  Z  9999998=0*  )  THEN  TUNC  HER:  DO; 

CALL  DENFMCH (FORMULA, A) ; 

PUT  FILE  (PUNCH)  LIST  (EOhKULA)  SKIP  ; 

END ; 

END  CARDS; 

CLOSE  FILE (PUNCH)  ; 

END  DISK; 

FORM  DETERMINANT : PROCEDURE; 

DO  1=1  TO  MODES  ; CET  (I)  ;  LET  (FOR  (I)  =M  (I,  I) -MU)  ; 

DO  .7  =  1  TO  MODES;  CE1  (J)  ; 

IF  J<I  THEN 

LET  (ROW  (I)  =  (M  (1,0)  ,10*1  (I)  )  )  ; 

IF  J>I  THEN  LET  (ROW  (i)  =  (i  OW  (I)  ,M  (I,  J)  )  )  ; 

END; 

END; 

END  FORM  DETERMINANT; 


DCL 

FORMULA  CHAR  (800) VAR, 

ITERATE#  FIXED  BINARY, 

JJOTEST  BI rJ  (1 )  INITIAL  ( 1  0  ’  B)  , 

MODES  FIXED  BINARY, 

DENFI1C3  EK'ii.Y  (FIXED  BIN  (31 , 0)  , FIXED  DIN  (31 ,0)  )  EXTERNAL, 
VALUE  CHAR  (IOC)  VA’}, 

PARK  PIXEL  BINARY, 

NOTHING  FIXED; 


END  OF  PROGRAM: 
END  FIRST; 


(CHECK  (EIGENVECTORS,  SOLVE, SEGREGAI  F.,  LEFI_  EIGEN  VEC'x  ORS)  )  : 
SECOND:  PROCEDURE  OPTIONS  ( .'1  A I  N )  ;  FOR  MAC_  OPTIONS  ; 


/*  JUNE  27,1  972  DIAGO  NALIZATIC  N  CE  LINEARIZED  COUFLED  MODE  EQUATIONS 
ADAPTtD  FROM  FORTIONS  OF  THE  CARLMAN  PROGRAM  */ 


/*  1 . */ 

/*2 . */ 

/*2A<MU*###*#s*1mm|»#V 
/*2B.  N  N  N  N  N  N  N  N  N  ti  N  N  N  N  N  *  / 

/*5 . */ 


CALI  INPUT; 

CALL  DETERMINANT (DIMENSION)  ; 
CALL  SEGREGATE; 

CALL  N UN ERIC AL_E VALUATION ; 
CALL  R  TGilT_  EIGEN  V ECT C R S  ; 

CALL  LEFT_ EIGEN VECTORS; 

CALL  EIGENVALUES; 


INPUT: PROCEDURE;  CN  ERCFILE  (S  YSIK)  GO  TO  EOF; 

GE'J  DATA  (NROW)  ;  DIKLNFIOfJ,  <SDI  N=KROH; 

DCL  PAIR  CHAR  (200) VAR ; 

GET  FILL  (PAIRS)  LIST  (PAIR)  COPY  ;  CRT  (PAIR);  CLOSE  FILE(PAIRS)* 
LET  (cOTIN  ES  =  0)  ; 

CL  "’  (aniH;NROW)  ;  LEI  (Q$  ( 0)  =  0;  T Y F E=  C)  ; 

LET  (PROD  (1)  =1)  ; 

let  (QS  (-  1)  - u  ;  TYPE--  1)  ; 

N  DI  M ,  a  D I K ,  N  R  0  W  =  D 1 M  E  M  S 1 0  N  ; 

C  F  T  ( 5)  E I  N  ;NDIK)  ; 

LET  (I  (d)DlK)  =1)  ; 

DO  J=1  TO  a)DIM;  HUMBER  ( 1 , J)  - J  ;  END; 

DO  1=0  TO  N ROW  ; C E  i.  ( I)  ;  LET  (  N  ( I )  =  0)  ;  EH D; 

CLOSE  FILE  (S Y SIN)  ;  OPLN  FILE  (SYSIN) TITLE ( • ROWS •)  ; 

DO  i= 1  ”0  NKOWjCEl  (I)  ;  GET  LIST  (RESULTS) COPY ;  FORM  (RESULTS) • 

DO  J=1  TO  N  H  0  W  ;  CET(J);  LET  (TERM  (I,J)=ARG  ( J  ,  R  O  W  ( I )  )  )  ; 

LET  (TERM  (I,  J)  =  E  V  A  L  (TERM  (I,  J)  ,  PAIR)  )  ; 

END;  END; 

PUT  FILL (SYSCP) EDIT  ( »  SECULAR  DETERMINANT •)  (SKIP, A)  ; 

GO  TO  EOF; 

OUTPUT  SENTRY;  CHAREX  (RL.SULTS=X  (TYPE,  Q$  (TYPE)  )  )  ; 

GO  TO  PUNCHER; 

OUTPUT  ft  :  ENTRY  (  NAME)  ;  DCL  NAME  CHAR  (*)  ;  CALL  DRNFMCH  (RESULTS, NAME)  • 
PUNCHER:  1  ’ 

PUT  FILE  (1  YSCP)  L±ST  (RESULT  f)  SKIP; 

RETURN ; 

PUT  FILE  (SYSCP)  EDIT  (RESULTS)  (SKIP,X(1),A)  ; 

EOF:  RETURN;  END  INPUT; 


DETERMINANT: PROC  FLUKE  (NDIM)  RECURSIVE; 
DECLARE  (NDIM,  I)  FIXED  BINARY ; 


LET  (fflUKBS-olTIfttS+l)  ; 

DCL  STIFFS  FIXED  EINARY;  2TIMES  =  I  NT  EGE  R  ( <j)TI  N  E3 )  ; 

ONE_D ; 

IF  ND J M= 1  THEN  DO; 

FT  (NUHFJER  =  "HUMDEK  (oiDIM,  1)  »)  ; 

fy- 


♦  i 
+  2 

♦  3 

♦  4 

♦  5 

♦  6 

♦  7 
+  8 
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+  13  ’ 

+  14 

+  15 
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+  17 

t  13 
f  19 
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«-  24 

«■  25 
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*-  23 
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♦  30 

+  3  1 

♦  32 

♦  33  f 

+  34 

♦  35 

«■  36 

♦  37 

+  38 

♦  39 

+  40 

+  41 

+  42 

♦  43 

+  4  4 

+  45 

+  46 

♦  47 

♦  48 

♦  49 

♦  50 

+  51 

♦  52 

+  53 

+  54 

♦  55 

♦  56 

♦  57 

♦  53 

♦  59 

♦  60 
♦  61 
+  62 


then  ocio  nLU 

CALL  OUTPUT; 

NULL:  LETraTIfl£S=aTIHEs-1) ; 

RETURN; 

END; 


LET^I^TIMES)^)1-10  (  INTEGER  <^DI»)  “iNTEGER  OTIttES)  ♦  1)  ;  LET(I=»I«); 

LET ( N=" NUMBER (STIMES, I) " ) ; 

LET  (PROD  (alTIMES* 1 )  =  PROD  (3TIMES)  ♦  TERfl  (STIHES ,N)  )  • 

IP  IDENT (TROD (STI MESt 1) ; 0)  THEN  GOTO  END_OF_CALCULATION ; 

L-0 ;  DO  J=1  TO  NDIfl;  IF  J-»=l  THEN  DO; 

L=L+ 1 ; NUMBER (<ST IMES  ♦  1 , L) =NUHBER (3TIMES, J)  ; END; END; 

CALL  DETERMINANT (A DIM-1) ; 

LET  (1  =  1  (ST IMES)  )  ; 

END_OF_CALCULATION: 

END  CALCULATE; 


LET  (d)TlMES  =  STIMES-1)  ; 
RETURN;  END  DETERMINANT; 


SIGN: PROC  RETURNS (CHAR  (4) )  ;  LET  (SIGN=*1)  • 

DCL  N$  (7)  FIXED  EINARY; 

DO  L#  =  1  TO  NROWjCLT  (LI)  ;  N#  =  INTEGER  (I  (LI)  )  ;  N$(LS)=NUKBER(L#/N#)  ;END; 
PUT  EDIT  (  (N$  (Ki)  LO  K$  =  1  TC  NROW)  )  (  10  F(3r0)  ); 

DO  L#*1  TO  NROW;  IF  H$(Lf)-,=L#  THEN  DO  J#=L#*1  TO  NROW 
IF  N5  ( J I) =L#  THEN  DO;  * 

N$(J«)  =Ni(L#)  ;  LET  (SIGN  =  -SIGN)  ;  END; 

END; 

END; 

RETURN ('SIGN')  ;  END  SIGN; 

EIGENVECTORS: PROCEDURE ; 

PIGHT_EIG EN VECTORS: ENTRY; 

PUV  FILE  (S YSCP) EDIT ( •  RIGHT  EIGENVECTORS')  (SKIP  k\  * 

EIGENVECTORS#: ENTRY ;  ' 

CET (LEFT) ; 

LET (X00=1) ; 

DO  J =2  TO  DIMENSION;  CET  (J)  ; 

7^.(D(J"1)=’X00*TtRM(J,1)  >  ;  USUALLY  X00-->1;  BAY  -->0  */ 

KAKF._5HALLER:  DO  1=2  TO  DIH;CET  (I)  ;  DO  J=2  TO  DIM-  CET  111  • 

m:  EI,D  "UE-S"*L“R: 

3DIM, N ROW = DIMENSION; 

CST(3CIM);  LET (I  (SDIM) =1)  ; 

IF  ON  THEN  DO; 

ON  =*»ON ; 

-JUT  EILE(SYSCP)EDIT('  DETEOBINANT  FOR  EIGENVECTORS*  w<;ktP  i\ 
LET (TYPE-0)  ;  CALL  DETERMINANT (NROW)  ;  END-  EIGENVECTORS •)  (SKIP, A) 

CALL  SOLVE; 
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♦  76 

♦  77 
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♦  103 
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♦  124 

♦  125 
;♦  125 

♦  127 

♦  128 


*  HOW,  3  cm,  DIMENSIGN=NR0W*1; 
end  EIGENVECTORS;  < 


LEFT  EIGENVECTORS rPROCE DU  RE; 

LEFT= 100 ; 

PUT  FILE  (SrSCP) EDIT ('  LEFT  EIGE N VECTORS ')  (SKIP, A)  : 
CET  (3DIM)  ; 

/*  FORM  THE  TRANSPOSED  MATRIX  */ 

DO  1=1  TO  NROW; CET (I) ; DO  J=1  TO  NFOW;CET(J); 
LET(TERH  (I,J)=EVAL  (AEG  (I,ROW(J)  )  ,PAIR))  ;  END;  END; 
CALL  EIGENVECTORS!; 

END  LE  FT_EIGEN VECTORS ; 


SOLVE: PROCEDURE (DIMENSION)  ;  /* 


INTERFACE  FOR 
FOR  SOLVING 


DO  J  =  1  TO  DIM;  CET  ( J) 
DO  1=1  TO  DIM: 


DCL  (1/ J)  FIXED  BINARY 
OPTS El (PRINT)  ; 

OPTSEI (NOPRINT)  ; 

V  ARIAELE_SOLUTICN 
REPLACE  A  COLUMN 
CET  (I)  ; 

LET  (TSAVE  (I)  =  TERH  (I,  J)  )  ; 

LET (TERM (I,J)  =D (I) )  ;  END  REPLACE 
LET  (TYPE=J«-LEFT)  ;  LET (Q$  (TYPE)=C) 
REST’ORE_PREVIOUS_  COLUMN: 

DO  1=1  TO  DIM  ; C£T  (I)  ;  LET  (TERM  (I,  J)  =TSAVE  (I)  ) 

end  variadle_sclution; 

EOSR-.END  SOLVE: 


CALLING  THE  DETERMINANT 
SIHUTLANEOUS  EQUATIONS*/ 


A_COLUMN  ; 

CALL  DETERMINANT  (3DIM) 


END 


OPTSET  (PRINT)  ; 


NUMERICAL_EVALUATION : PROCEDURE;  DCL  PAIR  CH A R (200) VA R • 
N  R  A  W=  N  ROW*- 1 ; 

EOF:  DO  1  =  0  TO  N  ROW  ;  CET(I);  F  $  ( I  ♦  1)  =  AR  ITH  (  P$  ( I)  ) 

PUT  LIST  ((PS (I)  DC  1=1  TO  NRA  W) )  ; 

LOSR:  RETURN;  END  NUMERICAL  EVALUATION; 


END: 


EIGENVALUES: PROCEDURE ; 

P$  =  P$/P$ ( N  R  A  W )  ; 

DO  1=1  TO  NROW;  Q» (I) =COMPLEX ( P  J  (NROW-I* 1 ) ,  0 . )  ;  END- 
DO  1=1  TO  NROW;  PUT  DATA  (Q# (I));  END; 

CALL  FRTC (Q*, NROW) ; 

IF  ER  ROR=  OK  THEN  PUT  LIST  ((Of  (I)  DO  1=1  TO  NPOW))SKIP; 

IF  ERROR-OK  THEN  PUT  FILE(SYSCP)  CATA((Q!(I)  DO  1=1  TO  NROW))SKIP- 
/*  FOB  EACH  EIGENVALUE*/  ’ 


DO  1=1  TO  NROW ;CET  (I)  ; 

FLOATA  (EIGEN (I) =  H E A L (Q#  (I) ) )  ;  F10ATA  (D=IMAG (Q# (I) ) ) 
LET(EIGFN(I)=EIGLN(I)  ♦!I*B)  ; 

LET  (MO  (I)  =  EIGEN  (I)  )  ;  CHARFX (STRlNG  =  HU  (J) )  ; 

CALL  OUTPUT!;  ATOMIZE  (MU (I) )  ; 

LET  (DET(I)=0)  ;  DO  J=1  TO  INTEGE  R (Q  $ (  0));  CET  (J)  ; 
LET  (DET  (I) =EVAL (X (  0,J)  RU, EIGE H ( I) ) »DET (I) ) ; 

BHD; 
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IF  IDEKT  (DEI  (I)  ; 0 )  THEn  DO;  ♦  195 

PUT  LIST ( •  EIGENVECTOR  ',1,  'SCI  FOUND');  ♦  195 

GOTO  EIGENHODES;  INC ;  ♦  197 

♦  198 

LEFT_RIGUT:  DO  L=0  to  1;  CET(L);  /•  L=0-->RIGHT,  L=1-->LEFT*/  ♦  199 

♦  200 

COHPo N ENTS :  DO  J=1  TO  NROW-1;  CET  (J)  ;  LET (JL=J*L*1 00)  ;  ♦  201 

LET  (CONDON  (I,  J,L)  =0)  ;  NTOTAL=INT£GEB (Q$ (JI )  );  ♦  202 

♦  203 

DO  K=  1  TO  NTOTAL;  CET(K);  ♦  204 

L£T(cCMPCN(I,J,L)=COHPCN(I,J,L)*EVAL(X(JL,K)  , HU, EIGER  (I)  ) )  ;  *205 

END,  _  -  ♦  206 

♦  207 

LET  (COUPON  (I,J,L)  =  CGflPON  (I,  J ,  L)  /DET(I)  )  ;  ♦  203 

LET  (COMP(I,J,L)=EVAL  (COHPON  (I,  J,L)  ,#I,-»I)  )  ;  ♦  209 

El  D  COHPO NENTS ;  ♦  210 

♦  211 
♦  212 

AGAIN:  DO  J=  NROW  TO  2  DY  *1;  CET(J);  LET(  ♦  213 

COUPON  (I,J,L)  =COMPOH  (I,  J-  1,L)  ;  CCHP  (I,  J  ,L)  =COKP  (I ,  J- 1 ,  L)  )  ;  END;  ♦  214 

L ET  (COKFON  (I,  1  f  L)  =1 ;  COMP  (I,  1  ,L)  =  1)  ;  '♦  215 

V  216 

LET(SC  =  0)  ;  ♦  217 

DO  J=1  TO  NROW-1;  CET  (J)  ;  LET  (SC=SQ*COMPOH  (I,J,I.)  *COHP  (I,  J,L) )  ;  END;  ♦  213 

Do  0=1  TO  NROW-1 ;CET (J)  ;LKT  (CCMFCN (I, J  #  L)  =COMPON (I, J,  L) /SQRT  (SQ) )  ;  END ; ♦  219 

♦  220 

END  LEFT_RIGHT;  ♦  221 

♦  222 

♦  223 

♦  224 

PUT  FILE  (SYSCP)  SKIP;  ♦  225 

DO  J=1  TO  NROW-1;  CET  (J)  ;  ♦  226 

LET  (S  (I ,  J)  -COUPON  (I , J  ,  0)  ;S1(I,J)=COMPON(I,J,1))  ;  ♦  227 

C II A  REX  (STRING=S  1  ( I,  J )  )  ;  PUT  FILE  (SYSCP)  LIST  (STRING)  ;  ♦  223 

CHAREX (STRING=S  (I,J));  PUT  FILE (SYSCp) LIST (SIR ING)  ;  ♦  229 

END;  ♦  230 

ElGEN  MODES :  ♦  231 

END;  ♦  232 

♦  233 

END  EIGENVALUES;  «.  234 

♦  235 

♦  236 

♦  237 

♦  238 

SKGREGA?E:PROCEDURE;  «.  239 

DO  1=0  TO  NROW* 1; CET  (I)  ;  LET (PS  (I ) =0) ;  END;  4  240 

NTERM5  =  I NTEGER (Q$ (-1 ) )  ;  DO  N$=1  TO  NTERMS;  CET  (NS)  ;  ♦  241 

IF  LOP  (X  (-  1 ,  NI) )  =  24  THEN  DO;  N  ST  A  RT=  1 ;  KSTOP=N»HGS  (X  (- 1#H$)  )  ;  END;  ♦  242" 

ELSE  NSTART,NSTOP=0;  .  ♦  243 

DO  N#=NSTART  TC  NSTOP;  CET  ( M)  ;  LET  (DU  M  P=A  RG  (N I  ,  X  (-  1 ,  N  $)  )  )  ; '  ♦'  244 

IF  IDENT (BUMF ;  0)  THEN  GOTC  HO  CONTRIBUTION;  \  ♦  245 

LET  (HP=H1GIIP0W  (  BUMP#  MU)  ;  PS  (HP)  =  P$  (IIP)  *EVAL  (BUMP,  HU,  1 )  )  ;  •  ♦  246 

NO_CONTRIBUTION:EKD;  *  247 

END  SEGREGATE;  «.  248 

♦  249 

♦  250 

DCL  DENFMC3  ENTRY  (FIXED  BI N  (3 1 , 0)  ,  FI  X ED  DIM  (3 1 r  0) )  EXTERNAL,  ♦  251 

(INPUT, OUTPUT)  ENTRY,  «.  252 

DIMENSION  FIXED  BINARY,  .  *  *  253 

(ALPHA, BETA)  FIXED  BINARY,  *  254 

SIGN  ENTRY  RETURNS  (CHAR  (4)  )  ,  .  «.  255 

RESULTS  CHAR (800)  VAR,  *  256 

NUMBER  (9,9) FIXED  El  NARY,  «.  257 

s) DIM  FIXED  BINARY,  «.  258 

NROW  FIXED  BINARY,  ♦  259 

DIM  FIXED  El N  DEFINED  DIMENSION,  ♦  260 

*» 

jr 


VALUE  CHAB  (100)  VAR, 

LEPT  PIXEL  EINARY <INII (0)  , 

nothing  fixed, 

P*(20) , 

DETER8INANT  ENTRY  (FIXED  BINARY)  • 

DCL  ON  BIT  (1)  I NI T  ( *  1  •  B)  ; 

DCL  STRING  CHAR(1C00)VAH; 

DCL  ERROR  CHAR  (1)  EXTERNAL, Q*  (20)  COHPLLX 
OK  CHAR  (1)  IN1T  ( • 0 • )  - 
END  SECOND; 


BINARY  PLOAT, 


♦  261 

♦  262 

♦  263 

♦  264 

♦  265 
,  ♦  266 
/  «■  267 

♦  268 

♦  269 

♦  270 


( 


»  SPART090. 7160Q 


//  EXEC  EXPORT, PROGPAS=EDITY  1  * ' T' K 6 5~G E? , T I D E- 
//,f,lKVL1B  DD  DSN=U*  P PL. LI  DRARY,  DISP=SHR 
//  UNIT=SYSDA^**^*^^  (THREE) , V0L=SER=RES103, 

//  DISP=  {OLD,  KEEP) 

//STSPRINT  DD  SYSOUT=A 
//SYS IN  DD  * 

•NEW  THIRD 

(CHECK (FORM_NEW_DRIVER, CHANGE  NON  LINEAR. 
L°prnc^0f<DL:R-S0I'UrI0N#  EVALUAI^-K0N  LINEAR  TERMS, 

FIHST_ORDER_ITeraiion, suastituteJintegrate 

ElGENVALUL’_CORRECIION)  )  :  n«ltl>KATK, 

THIRD: PROCEDURE  OPTIONS  (MAIN)  ;  PORMAC_OPTIONS ; 

/♦  NEW  VERSION  9-22-72  INCLUDES  FIRST  ORDER  ITERATION  ♦/ 


/*1*/  CALL  INPUT 
/*2*/  CALL  POR«_NEW_DRIVER; 


/*2. 5*/CALL  LOWEST  ORDER  SOLUTION* 
CALL  SAVER;  “  ~  ’ 

/♦3  */  CALL  FIRST_ORDER_ITER ATION* 
/*N*/  CALL  EV ALU AT E_NO H_LI NEAR  TERMS; 


CALL  CHANGE_NON_LINEAR; 


FORM_NEW_DRIVLR:  rifOcEDURE; 

DO  1  =  1  To  MODES ; CET  ( I)  ;  LET(SCK=C); 

DO  0=’  TO  MODES;CET(J)  ;LET(SUM=SUM*S(I,J)*D(J))  *END* 
LET  (NEWD{I)  =SUM)  .-ATOMIZE  (SUM)  •  M.END, 

END; 

DO  1=1  TO  MODES;  CET  (I)  ;ATOMIZE(D(I))  ;LET  (D  (I)  =  NEWD  (I) )  • 
ATOMIZE  (NEHD  (I)  )  ;  END*  *  M  * 

END  FORM  NEW  DRIVER;  ’  '  ’ 


LOWL’ST_ORDER_SOL  E TIO  N  :  PROCEDURE  * 

MODE: DO  1=1  TO  MCDES;  CET  (I); 

LP-LOP  (D  (I,  )  ;IF  L P  =  2 N  THEN  DO  ;NSTART=1;  NSTOP=NARGS  (D(I,  ,  ;FIJ 
LET (K$=0) ;  ”LS£  nSTART,NSTOP=0; LET (SUH=0 

TERMS: DO  N*  =  NSTART  TO  HSTOP; CET (Nt)  * 

LET  (TERB=AHG {H», D  (I) ) )  ; 

F  I  LENT (TERM; 0)  THEN  GO  TO  NODRIVE* 

LET(KXPOX=DKRIV  (LOG  (TERM)  ,TIHE)\’; 

/*  LET(TERH=TtRfl/(EXPOX*tI*KU.  (I)  )  ;  HOMOGEN=E  VAL  (TERM,!!  ME,  0) 

LET(K$  =  K$4-1  ;  LOW  (I,  K$)  =HCMOGEN*EXP(»I*MU.  (I)  *TIME)  )  ’  . 

LET  (K$=K$4l  ;  LOW  ( I ,  K$)  =110 HOGEN)  ;  k 

CALL  TRANSMOGRIFY { ’LOW (I,K$) •) • 

CHAREX (STRING*  LOW(I,K$))j  ATOMIZ E (LO W ( I  Kt\  \  • 

PUT  FILE  (HOMOGEN) LIST (•  * |  I STEZNG) SKIP*  '  '$)  1  * 

LF,T(B(0,I,K$)=LOW  (I,KS)  ;EXFONENT(0,’l,K$)=fI*MU  mi- 
SAVE(D(0,I,K$)  ;EXFONENT(0,I,K$))  •  '  '  fI 

ATOMIZE (TERM ; HOMOGEN) ; 

NOLRIVE: 

END  TERMS;  •  s 

LET(K$  (I)=K$)  ; 

ToTAL(0,I)=  INTEGER  (K$)  ; 

end  mode; 
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/•  CHANGE  9-26-72 

1  COT  DOWN  NUHEER  OP  TERNS  IN  THE  PINAL  ANSWER 

2  PFRMIT  LOW  SOLUTION  TO  APPEAR  EVEN  IF  DRIVER  IS  ZERO  IN 
DO  1=1  TO  NODES;  CET  ( I)  ;  LET ( K  $  1 ;  B  (0,  I  ,  1 )  =LOV  (I)  ;  LOW  (I  ,  1) 
K$  (I)  *1 ;)  ;  TOTAL  (0 /  I)  =1  ;  END; 


THAT  NODE*/ 
LOW  (I)  ; 


END  L ON EST_ORDE RESOLUTION ; 


EVALUATE^ NO N_LINEAfi_ TER MS: PROCEDURE; 
OPTSET (NOEDIT)  ; 

OPEN  FILE (PUNCH)  TITLE (' NLTERHS •) OUTPUT; 
OPTSET  (NOPRINT)  ; 


NODE :  DO  1=1  TO  NODES;  CET(I); 

LET  {I1A  HY  =  0)  ; 

DO  J=1  TO  NODES;  GET  (.1)  ; 

DO  K= 1  To  NODES ;  CET  (K)  ; 

IF  J  =  K  THEN  GOTO  NOSWEAT; 

IP  NOD=0  G  (  IDENT  (J  .-SPECIAL)  |  IEENT  (KjSPECIAI.)  )  THEN  GO  TO  NOSWEAT; 
NEST-J+K; 

NASTY=NEST/2  ;  IF  2*HASTY=NEST  THEN  GOTO  NOSWEAT; 
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IF  -.IDENT  (GANNA  (If  J,K)  ;  0)  THEN  CO; 

LET  (GAM=GANMA  (I,  J,  K)  )  ;  ATOMIZE  (GANliA  (I,  J,  K)  )  ; 

IF  LOP  (GAM)  =  24  THEN  DO;  NSTART--  1 ;  NSTOP=  NA  PGS  (G  AN) 
ELSE  NSTART, NSTOP=0  ; 

DO  NS=  NSTART  TO  NSTOP;  CET (NS); 

LET (GANMEH=ARG (NS,  GAM) )  ; 

DO  L$- 1  TO  INTEGER  (K$  (J) )  ;  CST(L$); 

DO  N$= 1  TO  INTEGER  (K$ (K) )  ;  CET (NS)  ; 

LET(MANY  =  !1ANY4  1  ; 

NONLINE (I, MANY) =GAMHER*LOW  (J,LS) *LOK (K, N$) )  ; 
PRIN T_OUT ( HO NLI N  t (I, MANY) )  ; 

CIIAREX  (STRING=NONI.lNF.  (I,  MANY)  )  ;  CALL  OUTPUT; 

SAVE (NONLINE (I, MANY) )  ; 

END;  END; 

ATOMIZE  (GAMMER)  ; 


END; 


END;  ATOMIZE (GAN) 


END; 

NOSWEAT; 
ENE  ;  END; 


LET  (  *#  (I)  =MAN  Y)  ; 
END  NODE; 


,.\V  " 


ITER=2;  CET (ITER) ;  LEt (K APPA  =  ITER ; K APP A 1  =  K APP A- 1 )  ;  KAPPA  1  =  ITER-V; 
PNINT_OUT (KAPPA; KAPPA  1)  ;  PUT  CAT A ( (TOTAL (2, IS)  DO  IS=1  TO  MODES)) 
CLOSE  FILE(PUNCH);  OPEN  FILE (PUNCH) TITLE (' Ti^EES ») OUTpUT ; 

INODES: DO  1=1  TO  NCDES ;  CET  (I)  ; 

DO  IB  =  1  TO  MODES;  CEI(IB);  LET  (1$  (ID)  =0)  ;  END; 

DO  J=1  TC  MODES;  CET (J) ; 


JMODE5 

KMODES 


THEN  GOTO 
TO  KAPPA  1 ; 


K_END; 

KSIG=KAPPA1 


DO  K  =  1  TC  MODES;  CET  ( K)  ; 

IF  IDENT  (GAMMA  (I, J,K) ;C) 

LOWE R_OR DEES :  DO  SIGHA-C 
CET (SIGMA ;KSIG) ; 

JTEFBS:  DO  JN=1  TO  TOTAL (SIGMA 
KTERMS;  DO  KN- 1  TO  TOTAL (KSIG 
OPTSET  (EXPND)  ; 

LET(EXPONA=  EXPONENT (SIGMA, . I, JH) ♦EXPONENT (KSIG,K,KN) )  ; 


SIGMA; 


,J) ; 

rK) 


CET  (JN)  ; 
CET  (KN) 
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Jp?si^0KL%EPL‘“|Eito"‘'"uus^  ; 

UR1SET (  PHINT) ; 

r«n^E=1  T0  n0D£S;  CET(IB);  ip  IDENT(EXPON;MU  (IB)  )  THEN  GOTO  NEW  DEAL 
E^D;  GOTO  NO_ACCOU  NT ;  NEW  DEAL-  t,UTU 

LET  (I$(lB)=Ii(lB)  *7; 

TEBH(I, IB,II (IB) ) =GAHHA  (I, J,K) *B(SIGnA,J, JN) *B  (KSIG,K,KN) 

*  E X F ( EX PON A) )  ; 


OP1SET (NOPRINT) ; 

ATOMIZE (EXPONA; EXPON)  ; 

N EW T  ER  M  I T  TH  T  c  ,  t  f ^ ^  ( * 1 G"' * ' " J " J >  :  EXPO  NENT  (  KSIG,  K  ,  !XN)  ; 
NE  TLR„(I,IB,i$(IBJ)  ;B(SIGMA,J,JN);B(KS1G,K,KN)  )  • 

CHARlX(STRING=NLWTERM(IflBfIi  (IB)))  ;  CALL  OUTPUT; 


no_account: 

END  KTEKMG; 

END  JTERMS ; 

END  LOWER  ORDERS; 

K_  EN  D : 

END  KHODES; 

END  J MODES ; 

NEW_TOTAL  (I)=1NTEgEr(I$)  ; 

END  IMODES; 

end  EVALUATE_NCN  linear  terms- 
OPTSET  (PRINT); 


EIG ENVALUE_CORR ECTION: PROCEDURE • 
OPTSET (ERSNP)  ; 

DO  1  =  1  TO  MODES;  CET  (I)  ; 

DO  J= 1  TO  INTEGER (#  t  (I)  )  ; 

CET(J)  ; 

OPTSET  (PRINT)  ; 

DO  J=1  TO  NEW  TOTAL  (I);  CET(J)- 
SAVE (NEWTERM  (I, J))  ; 

END; 

END; 

OPTS  ET (NOPRIN1)  ; 

END  EIGENVALUE_CORRECTION; 


TRANSMOGRIFY : PROCEDURE (A) ;  DCL  A  CHAR  (*)  - 
OPTSET  (  NOEX  PND)  ;  ' 

AT0h!z?U$$)  •  ET  (T0P=N(JM  {t$$)  1  WTTOB-DEHOfl  ($ J$,  )  ; 

PLSE°m<'T4  Rt°m  =  t  n  d  r,LN  D°:‘1S^BT=1;«STOP  =  NARGS  (BOTTOM);  END; 
ELSE  MSf ART, MSTOP=0 ;  LET ( NEWECT= 1) ;  1  * 

DO  M|  =  MSTART  TO  MSTOP ;CET  (M$)  ; 

LET (PACTOR  =  AHG (Mi,  BOTTOM) )  ; 

LET  (FAX  =  LVAL (FACTOR, #I,-*l) )  ; 

REALS;  LET  (NE WBOT  =EX P A N D  (  P ACTOR* P A  vi  * SiSWDOT  ;  T C  P=TO P* P A X )  ; 

OPTSET (PRINT) ; 

LEI  ("A'‘=TOP/NE*bCl)  ; 

ATOMIZE(NEWBOT;BOTlCH; TOP ;  FAXjFACIOR) ; 

EOS R ;  END  TRANSMOGRIFY ; 


OPTSET (NOPRINT)  ; 
CHANGE_NCN_LINEAR:  PROCEDURE; 

MODE; 

DO  1=1  TO  MODES;  CET  (I) ; 
DO  J=1  TO  MODES ;  CET  (J)  ; 
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CO  K= 1  TO  MODES ;  C£T  (K)  ; 
LET  (SUM=0)  ; 


DO  HU  =  1  TO  MODES;  CET(MU); 

DO  LAMBDA  =  1  TO  MODES;  CET(LAHBDA); 

DO  NU=  1  TO  MODES;  CET(NU); 

LET  (SUH=SUH  +  GAMMA  (MU, LAMBDA, NU)  *S1 (LAMBDA, J) *31  (NU,K)  * 
S (I, HU) )  ; 

END;  END;  END;  ATOMIZE(MU); 

LET (NEUG AM (I, J, K) =  SUM) ;  ATOMIZE (SUM) ; 

END  MODE; 


FANDANGO: 

DO  1=1  TO  MODES ; CET  (I)  ; DO  J=1  TC  MODES ; C ET  (J)  ; DO  K=1  TO  HODES;CE? (K)  ; 

LET (GAMMA (I,  J,K)  = NEW GAM  (I, J,K) )  ;  ATOMIZE ( N E WG A M ( I , J , K) )  ; 

END  FANDANGO; 

CALL  ATOMIC (' GAMMA (I, J,K) ' ,3) ; 

END  CHANGE_NCN_LINEAR ; 


DCL  ATOMIC  ENTRY (CHAR  (25) VAR, FIXED  BINARY); 

ATOMIC:  PROCEDURE  (A,  NE)  ;  DCL  A  CflAR(*),ND  FIXED  DIN; 

DCL  STRING  CHAR  (2000)  VAR, C  Ci!AR(20)  VAR,  IB  FIXED  BIN; 

TB=lNDEX  (A,  •  (')  -  1  ;  C-SUBSTR  (A,  1,  IB)  |  |  1  1  '  ;  OPEN  FI  LF  (SYSCP)  TITLE  (C) 
OUTPUT; 

DCL  NT  (4)  ;  NT  =  1 ; DC  NJ=1  TC  ND; NT (NJ) =HODES; END; 

ALL  :  DO  1=1  TO  NT  ( 1 )  ;  DO  >1  =  1  TO  N  T  ( 2 )  ;  CET(I;J); 

DO  K  =  1  TO  NT  (3)  ;  CET  (K)  ;  CO  L=  1  TO  NT  (4)  ;  C  E7  (L)  ; 

IF  -«D  ENFMCG  ( '  Z99  9  9  997=  •  |  |A,' 7.9999993=0*)  Then  DO; 

CALL  DENFMCil  (STRING,  A)  ;  CALL  DENFMC8  (A  |  |  '  =0  ' )  ; 

PUT  FILE (SYSCP) LIST (STRING) SNIP; 

END; 

END  ALL; 

CLOSE  FILE (SYSCP); 

RETURN; 

END  ATOMIC; 


INPUT: PROCEDURE; 

/*1.*/  GET  FILE (CONTROL) DATA  (MODES, 

SECOND_ORDER, 

NOTHING)  ; 

MOPAL=MODES/2 ;  IF  2* MODA L= MOD ES  THEN  NOD=1;  ELSE  NOD=0* 
TOTAL  =  0 ;  NEW  TOT  AL  =  0 ; 

OPTSET(EXPNU)’ ; 

/*2.V  /•  ZERO  THE  ARRAYS  OF  COEFFICIENTS  */ 

/•A*/  CALL  FCRMAC  ER AS E  ( • D (I) » , MOOES , 1) j 
/•B*/  CALL  EORMAC  EP ASt  ( • S 1 (I , J) • , MODFS, 2) ; 

/•CV  CALL  FuRMAC  ER  ASM  *  GAMMA  (I  ,  J  ,  K)  •  ,  MCOES ,  J)  • 
/•D*/  CALL  FORMAC_ERASl  ('S(I,J)  •  ,  MOL  hi. ,  ^  )  ; 

/•£•/  CALL  FORMAC_IRASl (• B(I,J)  • , MOLES, 2)  ; 

FORM AC_ ERASE: PROCEDURE (A, N, RANK) ; 

DCL  A  CHAR  (*),  (N, RANK,  NT  (4)  )  FIXED  BINARY; 
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DO  1=1  TO  14;  Ip  j.<  =  R  ANK  THEN  NT(I)  =  N;  ELSE  NT  ( I )  =  1  ;  end* 

ZEHO:  DO  1=1  1C  NT  (1)  ;CET  (1)  ;  CO  J=1  10  NT  (2 )  ; C^T 
DO  K-t  TO  NT(J);  CET  { K )  ;  DO  1=1  TO  N  T  ( 4 )  ;  CET(L); 

RE1  (A=0)  ;  END  ZE’SO;  ATOM  IZE  (I  ;  J  ;  K;  L)  ;  END  POPfiAC_ERASE; 

/* 3. */  /♦  INPUT  THE  ARRAYS*/ 

/♦A*/  CALL  EQUATIONS  (' DRIVER* )  ; 

/*B*/  call  equations  ('s •) ; 

/*C*/  CALL  EQUATIONS  ('GAMMAS')  ; 

EQUATIONS: PROCEDURE (A)  ;  DCL  A  CHAR/*)* 

OPEN  FILE  (SYSIN)  TITLE  (A)  INPUT; 

ON  ENDPILF (SYSIN)  ELGIN ;  CLOSE  FILE(SYSIN);  GOTO  EOF*  END- 

LOOP:  GuT  FILL (SYSIN)  LIST  (VALUE) COPY ;  FORM  (VALUE)  •  GOTO  LOOP- 
EOF:  RETURN;  END  EQUATIONS;  J  *  4U1U  LOOP, 

DO  I D=  1  TO  MODES  ;CET  (IB)  ;  L E 1  (B U  (I  D)  =E XP A N  D  ( -  «  I *MU  (I B )  )  )  •  END  - 
LET  (MULIS1’S  (MU.  (  1)  ,  MU  (1)  )  )  ;  '  '  ' '  ’ 

DO  IS  =  2  TO  M  DDES ;  CET  (IS); 

LET  (MULJ.ST=  (MULIST,  HU.  (IS)  ,  KU  (IS)  ));  end- 
PRI NT_OUT ( MU L 1ST ) ; 


/*<4.*/  /*  FORM  THE  SECOND  ORDER  EQUATIONS  FOR  RtFERENCF  ♦/ 

DCI  S  ECO N  D_CRDEH  B I T  <  1 )  INIT  ( »  1 '  1)  ; 

DCL  S  ECON  D_OR  CLH_  EQU  A 1  IONS  ENTRY;  IF  SECOND  ORDER  THEN  CAI 7 
SECOND  ORDER  EQUATIONS; 

ELSE  GOTO  LOSE; 

S ECCN D_0 R DER_ EQU Al ION S : PROC ;  DO  1=1  TO  MODES-  CET/I1- 
LET  (DRIVE  (I)  =  DERI  V  (C  (I)  ,  71 M  E)  )  ;  *  (  1  * 

DO  IS-1  TO  MODES;  CET  (IS)  ;  LET  (DRIVE  (I)  =DRIVE  (I)  +M  (I,  IS  )*D(  IS  )) 

■  w  5=-i  !o  ’em 

(I  )  .LET (LINtuR ( I,  J) =L1HEAR  (I ,  J) *M (I,  IS  )  *M (  IS  ,J))-END-  » 

j  j  ?0  K  =  1  CET(K);  LET(MATHIEU(I,J,K)=GAMMA  (1 , 0  ,  K)  *G  A  MM  A  ( I ,  K ,  J  ♦ 

LET (QUAD  (I, J, K) =0  );  CO  IS=1  TO  MODES;  CET  ( IS)  ; LET  (QUAD (I  J  K) -  i 

?S  ;e5d;(  “  'J'K,+"(  15  IS  ,  «>  ''♦«<♦ 

DO  L=  1  TO  MODES;  C  E  T  ( I )  ;  LF.T  (TR  IPLE  ( I ,  J,  K,  L)  =0  );  DO  IS=1  TO  Hum-* 

’j  lm^amm  w7  (I'J'K,L)  =iniPLf:  (1,‘1,K'L)  (I'  IS  »K)*GAHMM  IS  • 

•  J,L)  *GAMMA  (j.,  J,  IS)  *GAMMA  (IS,  K,  L)  )  ;  END;  1 

END;  END;  END;  4 

END  SECCND_ORDER_EquA1IONS;  * 


/♦5. */  /♦  PRINT  OUT  THE  TERMS  OF  SECOND  ORDER  */ 

CALL  FORMAC_ PRINT ( • DRIVE  (I)  •  ,  1)  ; 

CALL  FOR MAC.PRINl ('LINEAR  (I, J)  • ,2)  ; 

CALL  F0RMAC_PIUNT('QUAL>(I,J,K1  '  ,  3)  ’  ; 

CALL  F0RHAC_PKINT('MAXHIEU  (1,2, K)  *,3)  ; 

CALL  FORMAC_FRINT(' TRIPLE  (I,  J,  K,L)  •  ,4)"  ; 

r>u?"p;GE;INr:PB0C  <f,'HANK>  ;  tlCL  *  CIIAHC),  (SANK,  »T  («)  )  PIUD  HIM; 

NT=1;  DO  1=1  TO  RANK;  NT  (I) =MCDES l END- 
PRI NT_ ALL : 

DO  1  =  1  TO  NT  ( 1)  ;  CET  (I)  ;  DO  J  =  1  TO  NT  (2)  ;  CET  (J)  • 

DO  K  =  1  TO  NT  (3)  ;  CEI  (K)  ;  DO  L=  1  TO  NT  (4)  ;  CEt  (L)’ • 

IF  -'DEN  FMCG  ('29959997='  |  |  A,' Z  9999  99  3  =  0*  )  THEN 
CALL  DE  MFMC2 (A)  ;  CALL  DENFMC9  (A  |  | '=0 ')  ; 

ENL;  END;  END;  END  PRINT  ALL; 

END  PORHAC_PRINT  ; 

DCL  PORMAC_P3INT  ENT RY (CHAR  (20)  VAR, FIXED  BINARY) 

DCL  FORMAC_EHASE  ENTRY (CHAR  ( 15) VAR, PIXED  BIN, FIXED  BIN), 


* 

♦ 

« 

* 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

♦ 

* 

+ 

* 

f 

+ 

* 

f 

+ 

♦ 


253 

2514 

255 

25o 

257 

258 

259 

260 
2*1 
2b2 

26  3 
26<4 

265 

266 

267 

268 

269 

270 

27  1 

272 

273 

274 

27  5 

276 

277 

278 

279 

280 
281 
282 
20  3 
234 

28  5 
?3  6 
2-37 
2  3  8 
.  8  9 
270 
27  1 

2  72 
2  >  1 
27  4 
2'ir 
24;, 
?07 

29  8 
294 

300 

301 

302 

30  3 
304 

306 

3  0  ( 

307 

308 

309 
110 

311 

312 

313 

114 

-115 

316 

317 
316 


EQUATIONS  ENTRY  (CHAP  (6)  VAR)  ; 


GOTO  EOS R ;  OUTPUTiENTRY; 

PUT  FILE(PUNCH)LIST (STRING)  SKIP; 

EOSR:  END  INPUT; 


ITERATE:  PROCEDURE; 

FIRST_ORDLR_ITERATlCN  :  ENTRY  ;  IT  E RATEI=1 ; 

ITERATIONS:  DO  ITLh=  1  TO  ITERATE!;  CET (KAPPA=ITER)  ;  CET  (ITER)  ; 
LET  (KAPPA  1  =  KAPPA- 1)  ;  KAPPA 1  =  ITER-  1 ; 

MODE:  DO  1=1  TO  MODES;  CET(I); 

/♦I.*/  CALL  SUBSTITUTE; 

/*2. ♦/  CALL  INTEGRATE; 

END  MCDE; 

END  ITERATIONS; 

EOSR:  RETURN;  END  1TERATF ; 


SUBSTITUTE: PROCEDURE; 

LET  (14=0)  ; 

JMODES:  DO  J=1  TO  MODES;  CET  (J)  ; 

KMODES :  DO  K=  1  1C  MODES;  CrT(K); 

IF  I  DENT  (GAMMA (I, J, K)  ; 0)  THEN  GOTO  K_  END ; 

I  0WLR_0HDERS:  DO  S 1 G M  A  =  0  TO  KAPPA  1;  KSIG=KAPPA 1-SI3JU  ; 

C £1  (SIGMA  ; KSIG)  ; 

JTERMS:  DO  J  N  = 1  TO  TOTAL (SIGMA  ,J);  CET (ON)  ; 

KTERMS:  DO  KH=1  TO  TOTAL (KSIG  , K)  ;  CFT (K  N)  • 

LET ( 14= 14 ♦ 1 ; 

NEW  TERM  (14)  =  GAMMA  (I,  J,  K)  *B  (SIGMA,  J,JN)  *  B  (KSIG,  K  ,  KN)  ; 
EXPONENT (i$)  =EX PC NL NT  (SIGMA, J,. IN) ♦EXPONENT (KSIG, K, KN)  ; 


SAVE  (  EXPONENT  (1  4)  ;  EXPONENT  (SIGMA,  J,  JN)  ;  EX  PON  E  NT  (  KSI  G ,  K ,  KN)  ; 
NEKTERM  (IS)  ;B  (SIGMA,  J,JN)  ;I3  (KSIG, K,  KN)  ); 

end  ktermo; 

END  JTERMS; 

END  LOWER  ORDERS; 

K_END : 

END  KMODES; 

END  JMODES; 

NE’K_TOTAL  (I)  =INTEGER  (14)  ; 

EOSR:  RETURN;  END  SUBSTITUTE; 


INTEGRATE: PROCEDURE; 

LKt(EIOEH  =  MU.  (I)  **I;  J4=-  1 )  ; 

IKOClSiDO  IS=1  TO  NEW  TOTAL  (I);  CET  (IS); 

LET  (EXPCNA  =  EXPONENT (14)  -EIGEN)  ; 

IF  -IDF.NT  (LXPCNA  ;0)  THEN  DO; 

LFT ( J t  =  J $♦ 2  ; 

EXPONENT (ITER,  I, JS) -EXPONENT (1$)  ; 

EXPONENT  (ITER,  I,J$+1)  =  E 1 G  E  N  ; 

E  (ITER, I,  J$)  =NLr,TKHH  (IS)  /EUONA;  E  (ITER,  I,  JSf  1) 
IF  I  DENT (B (IT  Eh,  I ,  J$)  ; 0)  THEN  LET ( J$=J$-2)  ; 
ATOMIZE  (NEWTERM (14)  J EX  PON A ; T E l F C W  (1$)  ; EX PONENT ( I $)  ; 
END; 


=  -B  (ITER,I,JX)  ); 
)  ; 


♦  310 

♦  320 

♦  321 

♦  322 

♦  323 

♦  324 

♦  325 

♦  326 

♦  327 

♦  328 

♦  329 

♦  330 

♦  331 

♦  332 

♦  333 

♦  339 

♦  335 

♦  336 

♦  337 

♦  330 

♦  339 

♦  390 

♦  39  1 

♦  39  2 

♦  39  3 

♦  399 

♦  395 

♦  396 

♦  397 

♦  34  ft 

♦  39  9 

♦  350 

♦  351 

♦  352 

♦  353 

♦  354 

♦  365 

♦  356 

♦  357 

♦  358 

♦  359 

♦  360 

♦  361 

♦  362 

♦  363 

♦  364 

♦  36  5 

♦  366 

♦  36  7 

♦  36  H 

♦  369 

♦  370 

♦  37  1 

♦  37  2 

♦  373 

♦  374 

♦  37  5 

♦  376 

♦  377 

♦  .3  7  0 

♦  370 

♦  390 

♦  39  1 

♦  39  2 

♦  39  3 

♦  394 


ELSE  PUT  t  jbtivi'  • 

END  moots; 

TOTAL  (ITER,  I)  =  I  N  TLGE  R  ( J  I )  ♦  1  ; 

DO  Ilr!  T0  TOTAL  (ITER, I)  ;  CET(II);  SAVE  ( 

D(ITER,I,Ii)  ;  TE t  f  0  W (ITER, I, 1 1)  ;  EXPONENT (ITER, I, If)  )  ; 
EOSR:  RETURN ; 

END  INTEGRATE; 


SAVER;  P30C ;  DO  ORDL(i  =  0  TO  1;  CET  (ORDER)  ; 

DO  J=1  TO  MODES;  CLT(J)  ; 

DO  K  J  =  1  TO  TOTAL  (J, ORDER)  ; 

CFT  ( K$)  ; 

PRINTOUT  (B  (Order,  J,  K$)  )  ; 

CHAR  f  X  (STRING- I!  (ORDER,  J,  K$)  )  ;  HIT  FILE  (SAVED)  LIST  (STRING)  ;  END  SAVER; 

DCL  (FOH(1_NEW_DRIVER,  CHANGE  NCN  LI !l E AR  ,  1.0 W KST_ 0 R D LR_SO L UTI0N , 

FI  HSl_ORDKR_  ITER  ATION , SUBST  I  TOTE,  INTEGRATED  A  VF.R, 

E VALUATION  ON _LlNEAh_T ERRS, El G t  N V  ALU E_CORR ECT IC N)  ENTRY  ; 

DCL  TRANS KOG R 11  Y  ENTRY  (CHAR  (20)  VAR)  ; 

DECLARE 

VALUE  CHAR  (800)  VAR, 

NEW_TOTAL ( 100)  FIXED  BINARY, 

(SIGMA, ORDER)  FIXED  BIN, 

OUTI  UT  ENTRY, 

input  entry, 

TOTAL (R0,20)  FIXED  BIN, 

STRING  CHAR (10000)  VARYING, 

DENFKC3  ENTRY  (FIXED  E  I N  ( 3  1 , 0)  ,  FI  XI  D  BIN  (31,0)), 

NOTHING  FIXED; 

end_cf_proc;ham:  end  third; 

//  LX  EC  EXE ORT,  P R OGR A H= E DTT Y 
//STEPLIB  DD  D;>N=U.FRL.  LIBRARY,  IISP=GHR 
//SYS PRINT  DD  SYSOUT  =  A 

//IN  CD  D S H  =  U  .  ROSEN.  A R P A  (THREE)  ,  DISP=Sff  R,  UNIT=S YS I)A ,  VOL  =  S  ER  =  R ES  1  0  3 
//OU;  DD  S  YSOUT-A 
//SYSIN  DD  ♦ 


♦  385 

♦  336 

♦  387 

♦  383 

♦  33  8 

♦  330 

♦  391 

♦  392 

♦  393 

♦  394 

♦  39  5 

♦  396 

♦  39  7 

♦  39  3 

♦  399 

♦  400 

♦  401 

♦  402 
4  40  3 

♦  404 

4  4  06 
4  40t 
4  407 
4  4  03 

4  4  09 
4  4  10 
4  4  11 
4  4  12 
4  413 
4  4  14 

♦  41  r 
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